require(readstata13)
require(estimatr)
require(mlogit)
require(MASS)

#### prevalence of legacy politicians in postwar Japan ####
Reed.Smith.data <- read.dta13("Reed-Smith-JHRED-CANDIDATES.dta")

legacy.prop <- matrix(NA, 25, 6)
for (i in 1:25) {
  temp.data <- subset(Reed.Smith.data, yr == i & result > 0 & result < 4)
  legacy.prop[i, 1] <- mean(temp.data$pre_mp)
  legacy.prop[i, 2] <- mean(temp.data$pre_mp & 
                              (temp.data$predrelation == 1 | 
                                 temp.data$predrelation == 2 | 
                                 temp.data$predrelation == 6 | 
                                 temp.data$predrelation == 7))
  legacy.prop[i, 3] <- mean(temp.data$pre_mp & 
                              (temp.data$predrelation == 1 | 
                                 temp.data$predrelation == 6))
  legacy.prop[i, 4] <- mean(temp.data$seshu)
  legacy.prop[i, 5] <- mean(temp.data$seshu & 
                              (temp.data$predrelation == 1 | 
                                 temp.data$predrelation == 2 | 
                                 temp.data$predrelation == 6 | 
                                 temp.data$predrelation == 7))
  legacy.prop[i, 6] <- mean(temp.data$seshu & 
                              (temp.data$predrelation == 1 | 
                                 temp.data$predrelation == 6))
}
legacy.prop <- legacy.prop * 100

HOR.year.seq <- c(1947, 1949, 1952, 1953, 1955, 1958, 1960, 1963, 1967, 
                  1969, 1972, 1976, 1979, 1980, 1983, 1986, 1990, 1993, 
                  1996, 2000, 2003, 2005, 2009, 2012, 2014)

cairo_pdf("Figure_1.pdf", 
          width = 4.8, height = 2.5, pointsize = 7, family = "Helvetica")
par(mar = c(3.5, 3.5, 0.5, 0.5), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(1947, 2018), ylim = c(0, 35), 
     main = "", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
lines(HOR.year.seq, legacy.prop[, 1], type = "b", pch = 19)
lines(HOR.year.seq, legacy.prop[, 2], type = "b", pch = 1)
lines(HOR.year.seq, legacy.prop[, 3], type = "b", pch = 4)
lines(HOR.year.seq, legacy.prop[, 4], type = "b", pch = 19, col = "gray60")
lines(HOR.year.seq, legacy.prop[, 5], type = "b", pch = 1, col = "gray60")
lines(HOR.year.seq, legacy.prop[, 6], type = "b", pch = 4, col = "gray60")
text(2014, legacy.prop[25, 1], "Type 1", pos = 4)
text(2014, legacy.prop[25, 2], "Type 2", pos = 4)
text(2014, legacy.prop[25, 3], "Type 3", pos = 4)
text(2014, legacy.prop[25, 4], "Type 4", pos = 4)
text(2014, legacy.prop[25, 5] + 0.2, "Type 5", pos = 4)
text(2014, legacy.prop[25, 6] - 0.65, "Type 6", pos = 4)
axis(1, at = c(1947, seq(1950, 2010, 10), 2014), 
     labels = NA, lwd = 0.5)
mtext(c(1947, seq(50, 90, 10), "2000", 10, 14), side = 1, 
      at = c(1946, seq(1950, 2010, 10), 2014), line = 1)
mtext("Year", side = 1, at = 1980, line = 2.5)
mtext("Percentage", side = 2, line = 2.5)
axis(2, lwd = 0.5)
dev.off()

#### perceived prevalence and stereotypes data ####
## survey on legacy members
stereotype.data <- read.csv("replication_data_stereotype.csv")

# number of respondents
nrow(stereotype.data)

# age squared (for the multinomial logit model)
stereotype.data$age2 <- stereotype.data$age ^ 2
# middle education
stereotype.data$middle.edu <- (stereotype.data$education == 3) * 1
# higher education
stereotype.data$higher.edu <- ifelse(stereotype.data$education < 4, 0, 1)
# LDP support
stereotype.data$LDP.support <- (stereotype.data$partisanship == 1) * 1
# non-LDP right party support
stereotype.data$right.support <- (stereotype.data$partisanship %in% c(4, 6)) * 1
# left party support
stereotype.data$left.support <- (stereotype.data$partisanship %in% c(2, 3, 5, 8, 9)) * 1
# independent
stereotype.data$independent <- (stereotype.data$partisanship > 9) * 1
# ideological extremity
stereotype.data$extremity <- abs(stereotype.data$ideology - 3)
# dummy variable for having seen legacy members elected from their districts
stereotype.data$elected.yes <- (stereotype.data$elected == 1) * 1
# dummy variable for having not seen legacy members elected from their districts
stereotype.data$elected.no <- (stereotype.data$elected == 2) * 1
# deviation of each respondent's perceived percentage of legacy members from the actual value
stereotype.data$misperception <- abs(stereotype.data$prevalence - legacy.prop[25, 4])
# actual value of the percentage of legacy members in 2014
round(legacy.prop[25, 4], 1)

## survey on women members
load("replication_data_women.Rdata")
# number of respondents
length(women.perceived.prevalence)

#### perceived prevalence of legacy (and women) members ####
# mean perceived prevalence
round(mean(stereotype.data$prevalence))
# median perceived prevalence
median(stereotype.data$prevalence)

cairo_pdf("Figure_2.pdf", 
          width = 4.8, height = 2.5, pointsize = 7, family = "Helvetica")
layout(matrix(1:2, 1, 2))
par(mar = c(3.5, 3.5, 2.5, 1), lwd = 0.5)
hist(stereotype.data$prevalence, 
     breaks = 20, freq = FALSE, ylim = c(0, 0.05), 
     main = "Legacy members", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Estimated percentage", side = 1, line = 2.5)
mtext("Density", side = 2, line = 2.5)
par(mar = c(3.5, 3.5, 2.5, 1), lwd = 0.5)
hist(women.perceived.prevalence, 
     breaks = 20, freq = FALSE, ylim = c(0, 0.05), 
     main = "Women members", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Estimated percentage", side = 1, line = 2.5)
mtext("Density", side = 2, line = 2.5)
dev.off()

# whether respondents had ever had a legacy member of the HoR elected from their district
round(prop.table(table(stereotype.data$elected)), 2)

#### who perceived the presence of legacy politicians as higher or lower? ####
## Table 1
# Model 1
prevalence.lm.result.1 <- lm_robust(prevalence ~ gender + age + I(age ^ 2 / 10) + 
                                      middle.edu + higher.edu, 
                                    stereotype.data)
round(summary(prevalence.lm.result.1)$coefficients, 2)

# Model 2
prevalence.lm.result.2 <- lm_robust(prevalence ~ gender + age + I(age ^ 2 / 10) + 
                                      middle.edu + higher.edu + 
                                      LDP.support + right.support + left.support + 
                                      ideology + extremity, 
                                    stereotype.data)
round(summary(prevalence.lm.result.2)$coefficients, 2)

# Model 3
prevalence.lm.result.3 <- lm_robust(prevalence ~ gender + age + I(age ^ 2 / 10) + 
                                      middle.edu + higher.edu + 
                                      LDP.support + right.support + left.support + 
                                      ideology + extremity +
                                      trust + ex.efficacy + in.efficacy, 
                                    stereotype.data)
round(summary(prevalence.lm.result.3)$coefficients, 2)

# Model 4
prevalence.lm.result.4 <- lm_robust(prevalence ~ gender + age + I(age ^ 2 / 10) + 
                                      middle.edu + higher.edu + 
                                      LDP.support + right.support + left.support + 
                                      ideology + extremity +
                                      trust + ex.efficacy + in.efficacy + 
                                      elected.yes + elected.no, 
                                    stereotype.data)
round(summary(prevalence.lm.result.4)$coefficients, 2)

# Model 2 with a different baseline for partisanship (note 13)
prevalence.lm.result.2.2 <- lm_robust(prevalence ~ gender + age + I(age ^ 2 / 10) + 
                                        middle.edu + higher.edu + 
                                        right.support + left.support + independent + 
                                        ideology + extremity, 
                                      stereotype.data)
round(summary(prevalence.lm.result.2.2)$coefficients, 2)

# Model 2 with a different baseline for partisanship and without ideology and extremity (note 13)
prevalence.lm.result.2.3 <- lm_robust(prevalence ~ gender + age + I(age ^ 2 / 10) + 
                                        middle.edu + higher.edu + 
                                        right.support + left.support + independent, 
                                      stereotype.data)
round(summary(prevalence.lm.result.2.3)$coefficients, 2)

## expected perceived share of dynastic members by respondents' age (Appendix B.2)
# most frequent value of gender and education
table(stereotype.data$gender)
table(stereotype.data$education)
# prediction
new.data.age <- data.frame(gender = 1, age = 18:80, 
                           middle.edu = 0, higher.edu = 1)
predicted.values.prevalence.age <- predict(prevalence.lm.result.1, 
                                           new.data.age, 
                                           interval = "confidence")

cairo_pdf("Figure_A1.pdf", 
          width = 4, height = 3, pointsize = 8, family = "Helvetica")
par(mar = c(3.5, 3.5, 1, 1), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(18, 80), ylim = c(35, 55), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
lines(18:80, predicted.values.prevalence.age$fit[, 1])
lines(18:80, predicted.values.prevalence.age$fit[, 2], lty = 3)
lines(18:80, predicted.values.prevalence.age$fit[, 3], lty = 3)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
mtext("age", side = 1, line = 2.5)
mtext("Estimated percentage of legacy members", side = 2, line = 2.5)
dev.off()

## expected perceived share of dynastic members by respondents' ideological self-placement (Appendix B.3)
# mean age
mean(stereotype.data$age)
# prediction
new.data.ideology <- data.frame(gender = 1, age = 49, 
                                middle.edu = 0, higher.edu = 1, 
                                LDP.support = 0, right.support = 0, left.support = 0, 
                                ideology = 1:5, extremity = c(2, 1, 0, 1, 2))
predicted.values.prevalence.ideology <- predict(prevalence.lm.result.2, 
                                                new.data.ideology, 
                                                interval = "confidence")

cairo_pdf("Figure_A2.pdf", 
          width = 4, height = 3, pointsize = 8, family = "Helvetica")
par(mar = c(3.5, 3.5, 1, 1), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.5, 5.5), ylim = c(45, 65), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
points(1:5, predicted.values.prevalence.ideology$fit[, 1], pch = 19)
segments(1:5, predicted.values.prevalence.ideology$fit[, 2], 
         1:5, predicted.values.prevalence.ideology$fit[, 3])
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Ideological self-placement (progressive/left\u2013conservative/right)", side = 1, line = 2.5)
mtext("Estimated percentage of legacy members", side = 2, line = 2.5)
dev.off()

## analysis of misperceptions of the prevalence of legacy members (Appendix C)
cairo_pdf("Figure_A3.pdf", 
          width = 2.4, height = 2.5, pointsize = 7, family = "Helvetica")
par(mar = c(4.5, 3.5, 1, 1), lwd = 0.5)
hist(stereotype.data$misperception, 
     breaks = 20, freq = FALSE, xlim = c(0, 100), ylim = c(0, 0.03), 
     main = "", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Deviation of estimated percentage\nfrom the actual value", side = 1, line = 3.5)
mtext("Density", side = 2, line = 2.5)
dev.off()

# reanalysis of Table 1 (Table A.1)
misperception.lm.result.1 <- lm_robust(misperception ~ gender + age + I(age ^ 2 / 10) + 
                                         middle.edu + higher.edu, 
                                       stereotype.data)
round(summary(misperception.lm.result.1)$coefficients, 2)

misperception.lm.result.2 <- lm_robust(misperception ~ gender + age + I(age ^ 2 / 10) + 
                                         middle.edu + higher.edu + 
                                         LDP.support + right.support + left.support + 
                                         ideology + extremity, 
                                       stereotype.data)
round(summary(misperception.lm.result.2)$coefficients, 2)

misperception.lm.result.3 <- lm_robust(misperception ~ gender + age + I(age ^ 2 / 10) + 
                                         middle.edu + higher.edu + 
                                         LDP.support + right.support + left.support + 
                                         ideology + extremity +
                                         trust + ex.efficacy + in.efficacy, 
                                       stereotype.data)
round(summary(misperception.lm.result.3)$coefficients, 2)

misperception.lm.result.4 <- lm_robust(misperception ~ gender + age + I(age ^ 2 / 10) + 
                                         middle.edu + higher.edu + 
                                         LDP.support + right.support + left.support + 
                                         ideology + extremity +
                                         trust + ex.efficacy + in.efficacy + 
                                         elected.yes + elected.no, 
                                       stereotype.data)
round(summary(misperception.lm.result.4)$coefficients, 2)

misperception.lm.result.2.2 <- lm_robust(misperception ~ gender + age + I(age ^ 2 / 10) + 
                                           middle.edu + higher.edu + 
                                           right.support + left.support + independent + 
                                           ideology + extremity, 
                                         stereotype.data)
round(summary(misperception.lm.result.2.2)$coefficients, 2)

misperception.lm.result.2.3 <- lm_robust(misperception ~ gender + age + I(age ^ 2 / 10) + 
                                           middle.edu + higher.edu + 
                                           right.support + left.support + independent, 
                                         stereotype.data)
round(summary(misperception.lm.result.2.3)$coefficients, 2)

# expected perceived share of dynastic members by respondents' age
predicted.values.misperception.age <- predict(misperception.lm.result.1, 
                                              new.data.age, 
                                              interval = "confidence")

cairo_pdf("Figure_A4.pdf", 
          width = 4, height = 3, pointsize = 8, family = "Helvetica")
par(mar = c(3.5, 3.5, 1, 1), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(18, 80), ylim = c(20, 40), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
lines(18:80, predicted.values.misperception.age$fit[, 1])
lines(18:80, predicted.values.misperception.age$fit[, 2], lty = 3)
lines(18:80, predicted.values.misperception.age$fit[, 3], lty = 3)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
mtext("age", side = 1, line = 2.5)
mtext("Estimated percentage of legacy members", side = 2, line = 2.5)
dev.off()

# expected perceived share of dynastic members by respondents' ideological self-placement
predicted.values.misperception.ideology <- predict(misperception.lm.result.2, 
                                                   new.data.ideology, 
                                                   interval = "confidence")

cairo_pdf("Figure_A5.pdf", 
          width = 4, height = 3, pointsize = 8, family = "Helvetica")
par(mar = c(3.5, 3.5, 1, 1), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.5, 5.5), ylim = c(30, 50), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
points(1:5, predicted.values.misperception.ideology$fit[, 1], pch = 19)
segments(1:5, predicted.values.misperception.ideology$fit[, 2], 
         1:5, predicted.values.misperception.ideology$fit[, 3])
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Ideological self-placement (progressive/left\u2013conservative/right)", side = 1, line = 2.5)
mtext("Estimated percentage of legacy members", side = 2, line = 2.5)
dev.off()

# percentage of overestimation
round(mean(stereotype.data$prevalence - legacy.prop[25, 4] > 0), 3)

# correlation between misperception and the original outcome
round(cor(stereotype.data$misperception, stereotype.data$prevalence), 3)

#### trait stereotypes about legacy members ####
# function to compute standard errors of the difference in proportion
prop.se <- function(p1, p2, n) {
  sqrt((p1 * (1 - p1) + p2 * (1 - p2) - 2 * p1 * p2) / n)
}

trait.labels <- c("Trustworthy", "Decisive", 
                  "Competent", "Politically experienced", 
                  "Honest", "Has strong leadership", 
                  "Benefits his or her\nconstituency", 
                  "Does not commit\ncorruption", 
                  "Wealthy", "Highly educated", 
                  "Has a broad network in\npolitics and business", 
                  "Likely to become\na cabinet minister")

# cumulative proportion
trait.stereotypes <- apply(stereotype.data[, 2:13], 2, 
                           function(x) prop.table(table(x)))
colnames(trait.stereotypes) <- trait.labels

# difference between "more applicable to legacy members" and "more applicable to non-legacy members"
trait.differences <- trait.stereotypes[1, ] - trait.stereotypes[3, ]
trait.order <- order(trait.differences)
trait.diff.se <- prop.se(trait.stereotypes[1, ], trait.stereotypes[3, ], 
                         nrow(stereotype.data))

cairo_pdf("Figure_3.pdf", 
          width = 4.8, height = 4, pointsize = 7, family = "Helvetica")
par(mar = c(2.5, 10, 1.5, 1.2), lwd = 0.5)
trait.barplot <- barplot(trait.stereotypes[, trait.order], horiz = TRUE, xlim = c(-0.01, 2.2), 
                         xlab = "", ylab = "", xaxt = "n", las = 1)
text(trait.stereotypes[1, trait.order[12]] / 2, 
     trait.barplot[12], "More applicable to\nlegacy members", col = "white", cex = 0.6)
text(trait.stereotypes[1, trait.order[6]] + trait.stereotypes[2, trait.order[6]] / 2, 
     trait.barplot[6], "No difference", cex = 0.6)
text(sum(trait.stereotypes[1:2, trait.order[1]]) + trait.stereotypes[3, trait.order[1]] / 2, 
     trait.barplot[1], "More applicable\nto non-legacy\nmembers", cex = 0.6)
mtext("Cumulative proportion", at = 0.5, line = 0, font = 2, cex = 1.2)
axis(1, at = seq(0, 1, 0.2), lwd = 0.5)
segments(seq(1.2, 2.2, 1 / 6)[3], 0, seq(1.2, 2.2, 1 / 6)[3], 14.4, 
         col = "gray50")
segments(seq(1.2, 2.2, 1 / 6)[-3], 0, seq(1.2, 2.2, 1 / 6)[-3], 14.4, 
         col = "gray50", lty = 3)
points(trait.differences[trait.order] / 1.2 + 4.6 / 3, trait.barplot, pch = 20)
segments(trait.differences[trait.order] / 1.2 + 4.6 / 3 + trait.diff.se[trait.order] * qnorm(0.025) / 1.2, trait.barplot, 
         trait.differences[trait.order] / 1.2 + 4.6 / 3 + trait.diff.se[trait.order] * qnorm(0.975) / 1.2, trait.barplot)
mtext("Difference", at = 1.7, line = 0, font = 2, cex = 1.2)
axis(1, at = seq(1.2, 2.2, 1 / 6), 
     labels = sprintf("%3.1f", seq(-0.4, 0.8, 0.2)), lwd = 0.5)
dev.off()

#### issue stereotypes about legacy members ####
issue.labels <- c("Education", "Crime and public\nsecurity", 
                  "Medical care", "Child welfare", 
                  "National security", "Declining birthrate", 
                  "Fiscal deficit", "Diplomacy", 
                  "Industrial policy", "Public works")

# cumulative proportion
issue.stereotypes <- apply(stereotype.data[, 14:23], 2, function(x) prop.table(table(x)))
colnames(issue.stereotypes) <- issue.labels

# difference between "legacy members are better" and "non-legacy members are better"
issue.differences <- issue.stereotypes[1, ] - issue.stereotypes[3, ]
issue.order <- order(issue.differences)
issue.diff.se <- prop.se(issue.stereotypes[1, ], issue.stereotypes[3, ], 
                         nrow(stereotype.data))

cairo_pdf("Figure_4.pdf", 
          width = 4.8, height = 3.5, pointsize = 7, family = "Helvetica")
par(mar = c(2.5, 10, 1.5, 1.2), lwd = 0.5)
issue.barplot <- barplot(issue.stereotypes[, issue.order], horiz = TRUE, xlim = c(-0.01, 2.2), 
                         xlab = "", ylab = "", xaxt = "n", las = 1)
text(issue.stereotypes[1, issue.order[10]] / 2, 
     issue.barplot[10], "legacy\nmembers\nare better", col = "white", cex = 0.6)
text(issue.stereotypes[1, issue.order[5]] + issue.stereotypes[2, issue.order[5]] / 2, 
     issue.barplot[5], "No difference", cex = 0.6)
text(sum(issue.stereotypes[1:2, issue.order[1]]) + issue.stereotypes[3, issue.order[1]] / 2, 
     issue.barplot[1], "Non-legacy\nmembers are\nbetter", cex = 0.6)
mtext("Cumulative proportion", at = 0.5, line = 0, font = 2, cex = 1.2)
axis(1, at = seq(0, 1, 0.2), lwd = 0.5)
segments(seq(1.2, 2.2, 0.25)[3], 0, seq(1.2, 2.2, 0.25)[3], 12, 
         col = "gray50")
segments(seq(1.2, 2.2, 0.25)[-3], 0, seq(1.2, 2.2, 0.25)[-3], 12, 
         col = "gray50", lty = 3)
points(issue.differences[issue.order] * 2 + 1.7, issue.barplot, pch = 20)
segments(issue.differences[issue.order] * 2 + 1.7 + issue.diff.se[issue.order] * qnorm(0.025) * 2, issue.barplot, 
         issue.differences[issue.order] * 2 + 1.7 + issue.diff.se[issue.order] * qnorm(0.975) * 2, issue.barplot)
mtext("Difference", at = 1.7, line = 0, font = 2, cex = 1.2)
axis(1, at = seq(1.2, 2.2, 0.25), 
     labels = sprintf("%5.3f", seq(-0.25, 0.25, 0.125)), lwd = 0.5)
dev.off()

#### analysis of heterogeneity in stereotypes about legacy politicians ####
# function to post-estimation simulation for the multinomial logit model (MNL)
prob.sim <- function(sim.data, par.draw) {
  exp.2 <- exp(tcrossprod(sim.data, par.draw[, seq(1, 11, 2)]))
  exp.3 <- exp(tcrossprod(sim.data, par.draw[, seq(2, 12, 2)]))
  denom <- (1 + exp.2 + exp.3)
  prob.1 <- colMeans(1 / denom)
  prob.2 <- colMeans(exp.2 / denom)
  prob.3 <- colMeans(exp.3 / denom)
  cbind(c(mean(prob.1), mean(prob.2), mean(prob.3)), 
        c(quantile(prob.1, 0.025), quantile(prob.2, 0.025), quantile(prob.3, 0.025)), 
        c(quantile(prob.1, 0.975), quantile(prob.2, 0.975), quantile(prob.3, 0.975)))
}

# function to estimate the MNL
MNL <- function(y, seed) {
  set.seed(seed)
  MNL.data <- mlogit.data(data.frame(y = y, 
                                     stereotype.data[, c("gender", 
                                                         "age", 
                                                         "age2", 
                                                         "middle.edu", 
                                                         "higher.edu")]), 
                          choice = "y", shape = "wide")
  result <- mlogit(y ~ 1 | gender + age + age2 + middle.edu + higher.edu, 
                   MNL.data)
  par.draw <- mvrnorm(1000, result$coefficients, vcov(result))
  original.data <- cbind(1, 
                         as.matrix(stereotype.data[, c("gender", 
                                                       "age", 
                                                       "age2", 
                                                       "middle.edu", 
                                                       "higher.edu")]))
  out <- list()
  sim.data <- original.data
  sim.data[, "gender"] <- 0
  out$male <- prob.sim(sim.data, par.draw)
  sim.data <- original.data
  sim.data[, "gender"] <- 1
  out$female <- prob.sim(sim.data, par.draw)
  out$age <- array(NA, c(3, 3, 63))
  for (age in 18:80) {
    sim.data <- original.data
    sim.data[, "age"] <- age
    sim.data[, "age2"] <- age ^ 2
    out$age[, , age - 17] <- prob.sim(sim.data, par.draw)
  }
  sim.data <- original.data
  sim.data[, "middle.edu"] <- 0
  sim.data[, "higher.edu"] <- 0
  out$low.edu <- prob.sim(sim.data, par.draw)
  sim.data <- original.data
  sim.data[, "middle.edu"] <- 1
  sim.data[, "higher.edu"] <- 0
  out$middle.edu <- prob.sim(sim.data, par.draw)
  sim.data <- original.data
  sim.data[, "middle.edu"] <- 0
  sim.data[, "higher.edu"] <- 1
  out$higher.edu <- prob.sim(sim.data, par.draw)
  out
}

# MNL estimation for trait stereotypes
trait.MNL <- list()
for (i in 1:12) {
  trait.MNL[[i]] <- MNL(stereotype.data[, i + 1], 12345)
  cat("Finished at", date(), "\n")
}

# MNL estimation for issue stereotypes
issue.MNL <- list()
for (i in 1:10) {
  issue.MNL[[i]] <- MNL(stereotype.data[, i + 13], 12345)
  cat("Finished at", date(), "\n")
}

cairo_pdf("Figure_A6.pdf", 
          width = 6, height = 8, pointsize = 8, family = "Helvetica")
layout(matrix(c(1:23, 23), 6, 4, byrow = TRUE))
par(mar = c(2, 2, 2.5, 0.5), lwd = 0.5)
for (i in 12:1) {
  plot(NULL, NULL, type = "n", bty = "n", xlim = c(0, 3), ylim = c(0, 1), 
       main = trait.labels[trait.order[i]], xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(h = seq(0, 1, 0.2), lty = 3, col = "gray50")
  points(seq(0.8, 1.2, 0.2), trait.MNL[[trait.order[i]]]$male[, 1], pch = c(19, 4, 15))
  segments(seq(0.8, 1.2, 0.2), trait.MNL[[trait.order[i]]]$male[, 2], 
           seq(0.8, 1.2, 0.2), trait.MNL[[trait.order[i]]]$male[, 3])
  points(seq(1.8, 2.2, 0.2), trait.MNL[[trait.order[i]]]$female[, 1], pch = c(19, 4, 15))
  segments(seq(1.8, 2.2, 0.2), trait.MNL[[trait.order[i]]]$female[, 2], 
           seq(1.8, 2.2, 0.2), trait.MNL[[trait.order[i]]]$female[, 3])
  mtext(c("Male", "Female"), side = 1, at = 1:2, line = 0, cex = 0.7)
  axis(2, lwd = 0.5)
}
for (i in 10:1) {
  plot(NULL, NULL, type = "n", bty = "n", xlim = c(0, 3), ylim = c(0, 1), 
       main = issue.labels[issue.order[i]], xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(h = seq(0, 1, 0.2), lty = 3, col = "gray50")
  points(seq(0.8, 1.2, 0.2), issue.MNL[[issue.order[i]]]$male[, 1], pch = c(19, 4, 15))
  segments(seq(0.8, 1.2, 0.2), issue.MNL[[issue.order[i]]]$male[, 2], 
           seq(0.8, 1.2, 0.2), issue.MNL[[issue.order[i]]]$male[, 3])
  points(seq(1.8, 2.2, 0.2), issue.MNL[[issue.order[i]]]$female[, 1], pch = c(19, 4, 15))
  segments(seq(1.8, 2.2, 0.2), issue.MNL[[issue.order[i]]]$female[, 2], 
           seq(1.8, 2.2, 0.2), issue.MNL[[issue.order[i]]]$female[, 3])
  mtext(c("Male", "Female"), side = 1, at = 1:2, line = 0, cex = 0.7)
  axis(2, lwd = 0.5)
}
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0, 3), ylim = c(0, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
legend("center", 
       legend = c("More applicable to legacy members\nlegacy members are better", 
                  "No difference between legacy and non-legacy members", 
                  "More applicable to non-legacy members\nNon-legacy members are better"), 
       pch = c(19, 4, 15), bty = "n", cex = 1.2, y.intersp = 1.5)
dev.off()

cairo_pdf("Figure_A7.pdf", 
          width = 6, height = 8, pointsize = 8, family = "Helvetica")
layout(matrix(c(1:23, 23), 6, 4, byrow = TRUE))
par(mar = c(2, 2, 2.5, 0.5), lwd = 0.5)
for (i in 12:1) {
  plot(NULL, NULL, type = "n", bty = "n", xlim = c(18, 80), ylim = c(0, 1), 
       main = trait.labels[trait.order[i]], xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(h = seq(0, 1, 0.2), lty = 3, col = "gray50")
  polygon(c(18:80, 80:18), 
          c(trait.MNL[[trait.order[i]]]$age[1, 2, ], rev(trait.MNL[[trait.order[i]]]$age[1, 3, ])), 
          border = NA, col = "#B0B0B080")
  lines(18:80, trait.MNL[[trait.order[i]]]$age[1, 1, ])
  polygon(c(18:80, 80:18), 
          c(trait.MNL[[trait.order[i]]]$age[2, 2, ], rev(trait.MNL[[trait.order[i]]]$age[2, 3, ])), 
          border = NA, col = "#B0B0B080")
  lines(18:80, trait.MNL[[trait.order[i]]]$age[2, 1, ], lty = 3)
  polygon(c(18:80, 80:18), 
          c(trait.MNL[[trait.order[i]]]$age[3, 2, ], rev(trait.MNL[[trait.order[i]]]$age[3, 3, ])), 
          border = NA, col = "#B0B0B080")
  lines(18:80, trait.MNL[[trait.order[i]]]$age[3, 1, ], lty = 2)
  axis(1, lwd = 0.5)
  axis(2, lwd = 0.5)
}
for (i in 10:1) {
  plot(NULL, NULL, type = "n", bty = "n", xlim = c(18, 80), ylim = c(0, 1), 
       main = issue.labels[issue.order[i]], xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(h = seq(0, 1, 0.2), lty = 3, col = "gray50")
  polygon(c(18:80, 80:18), 
          c(issue.MNL[[issue.order[i]]]$age[1, 2, ], rev(issue.MNL[[issue.order[i]]]$age[1, 3, ])), 
          border = NA, col = "#B0B0B080")
  lines(18:80, issue.MNL[[issue.order[i]]]$age[1, 1, ])
  polygon(c(18:80, 80:18), 
          c(issue.MNL[[issue.order[i]]]$age[2, 2, ], rev(issue.MNL[[issue.order[i]]]$age[2, 3, ])), 
          border = NA, col = "#B0B0B080")
  lines(18:80, issue.MNL[[issue.order[i]]]$age[2, 1, ], lty = 3)
  polygon(c(18:80, 80:18), 
          c(issue.MNL[[issue.order[i]]]$age[3, 2, ], rev(issue.MNL[[issue.order[i]]]$age[3, 3, ])), 
          border = NA, col = "#B0B0B080")
  lines(18:80, issue.MNL[[issue.order[i]]]$age[3, 1, ], lty = 2)
  axis(1, lwd = 0.5)
  axis(2, lwd = 0.5)
}
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0, 3), ylim = c(0, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
legend("center", 
       legend = c("More applicable to legacy members\nlegacy members are better", 
                  "No difference between legacy and non-legacy members", 
                  "More applicable to non-legacy members\nNon-legacy members are better"), 
       lty = c(1, 3, 2), bty = "n", cex = 1.2, y.intersp = 1.5)
dev.off()

cairo_pdf("Figure_A8.pdf", 
          width = 6, height = 8, pointsize = 8, family = "Helvetica")
layout(matrix(c(1:23, 23), 6, 4, byrow = TRUE))
par(mar = c(2, 2, 2.5, 0.5), lwd = 0.5)
for (i in 12:1) {
  plot(NULL, NULL, type = "n", bty = "n", xlim = c(0, 4), ylim = c(0, 1), 
       main = trait.labels[trait.order[i]], xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(h = seq(0, 1, 0.2), lty = 3, col = "gray50")
  points(seq(0.8, 1.2, 0.2), trait.MNL[[trait.order[i]]]$low.edu[, 1], pch = c(19, 4, 15))
  segments(seq(0.8, 1.2, 0.2), trait.MNL[[trait.order[i]]]$low.edu[, 2], 
           seq(0.8, 1.2, 0.2), trait.MNL[[trait.order[i]]]$low.edu[, 3])
  points(seq(1.8, 2.2, 0.2), trait.MNL[[trait.order[i]]]$middle.edu[, 1], pch = c(19, 4, 15))
  segments(seq(1.8, 2.2, 0.2), trait.MNL[[trait.order[i]]]$middle.edu[, 2], 
           seq(1.8, 2.2, 0.2), trait.MNL[[trait.order[i]]]$middle.edu[, 3])
  points(seq(2.8, 3.2, 0.2), trait.MNL[[trait.order[i]]]$higher.edu[, 1], pch = c(19, 4, 15))
  segments(seq(2.8, 3.2, 0.2), trait.MNL[[trait.order[i]]]$higher.edu[, 2], 
           seq(2.8, 3.2, 0.2), trait.MNL[[trait.order[i]]]$higher.edu[, 3])
  mtext(c("Low", "Middle", "High"), side = 1, at = 1:3, line = 0, cex = 0.7)
  axis(2, lwd = 0.5)
}
for (i in 10:1) {
  plot(NULL, NULL, type = "n", bty = "n", xlim = c(0, 4), ylim = c(0, 1), 
       main = issue.labels[issue.order[i]], xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(h = seq(0, 1, 0.2), lty = 3, col = "gray50")
  points(seq(0.8, 1.2, 0.2), issue.MNL[[issue.order[i]]]$low.edu[, 1], pch = c(19, 4, 15))
  segments(seq(0.8, 1.2, 0.2), issue.MNL[[issue.order[i]]]$low.edu[, 2], 
           seq(0.8, 1.2, 0.2), issue.MNL[[issue.order[i]]]$low.edu[, 3])
  points(seq(1.8, 2.2, 0.2), issue.MNL[[issue.order[i]]]$middle.edu[, 1], pch = c(19, 4, 15))
  segments(seq(1.8, 2.2, 0.2), issue.MNL[[issue.order[i]]]$middle.edu[, 2], 
           seq(1.8, 2.2, 0.2), issue.MNL[[issue.order[i]]]$middle.edu[, 3])
  points(seq(2.8, 3.2, 0.2), issue.MNL[[issue.order[i]]]$higher.edu[, 1], pch = c(19, 4, 15))
  segments(seq(2.8, 3.2, 0.2), issue.MNL[[issue.order[i]]]$higher.edu[, 2], 
           seq(2.8, 3.2, 0.2), issue.MNL[[issue.order[i]]]$higher.edu[, 3])
  mtext(c("Low", "Middle", "High"), side = 1, at = 1:3, line = 0, cex = 0.7)
  axis(2, lwd = 0.5)
}
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0, 3), ylim = c(0, 1), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
legend("center", 
       legend = c("More applicable to legacy members\nlegacy members are better", 
                  "No difference between legacy and non-legacy members", 
                  "More applicable to non-legacy members\nNon-legacy members are better"), 
       pch = c(19, 4, 15), bty = "n", cex = 1.2, y.intersp = 1.5)
dev.off()